The Microscopic Mechanisms of Nonlinear Rectification on Si-MOSFETs Terahertz Detector

Studying the nonlinear photoresponse of different materials, including III-V semiconductors, two-dimensional materials and many others, is attracting burgeoning interest in the terahertz (THz) field. Especially, developing field-effect transistor (FET)-based THz detectors with preferred nonlinear plasma-wave mechanisms in terms of high sensitivity, compactness and low cost is a high priority for advancing performance imaging or communication systems in daily life. However, as THz detectors continue to shrink in size, the impact of the hot-electron effect on device performance is impossible to ignore, and the physical process of THz conversion remains elusive. To reveal the underlying microscopic mechanisms, we have implemented drift-diffusion/hydrodynamic models via a self-consistent finite-element solution to understand the dynamics of carriers at the channel and the device structure dependence. By considering the hot-electron effect and doping dependence in our model, the competitive behavior between the nonlinear rectification and hot electron-induced photothermoelectric effect is clearly presented, and it is found that the optimized source doping concentrations can be utilized to reduce the hot-electron effect on the devices. Our results not only provide guidance for further device optimization but can also be extended to other novel electronic systems for studying THz nonlinear rectification.


Introduction
The terahertz (THz) wave, ranging from 0.1 to 10 THz, is a region of the electromagnetic spectrum between infrared and millimeter-wave bands and hosts the advantages of optics and microwave electronics [1]. The rising interest in THz technology in fields such as radar, communication, non-destructive detection, biomedical and environmental detection [2,3] has prompted the progressive development of compact THz sources and high-performance detectors.
Recently, there have been a number of positive developments in international research around the optoelectronic properties of electronic systems in low-dimensional materials. As the most classical two-dimensional material, it has shown new vitality in a series of important applications such as ultra-broadband, high-speed optical communication and biosensing, which overturn traditional functions [28]. Meanwhile, other low-dimensional materials, such as graphene [29][30][31], black phosphorus [32][33][34] and transition metal dichalcogenides (TMDCs) [35][36][37][38][39], all exhibit extremely strong nonlinear THz rectification effects. Many efforts have been made to realize the modulation of the structural and electronic properties of low-dimensional semimetal materials by optical techniques, such as high-order harmonic generation [40] and light frequency mixing [41]. These works open new possibilities for the development of dissipation-less and ultrafast topological devices in data processing, sensing, and communication [42]. However, the performance of THz photodetectors is limited due to the low photon energy of the THz band, weak light absorption, and low collection efficiency of carriers [43]. In addition, conventional THz photodetectors are also limited by low speed, low operation temperature, and high power requirements. Therefore, it is challenging to realize room-temperature high-performance THz photodetectors [44,45]. To overcome these difficulties, we have explored the physical process of nonlinear THz rectification based on a typical material-silicon. It will be of great help for further research on future THz nonlinear rectification mechanisms based on graphene, black phosphorus, topological Dirac/Weyl semimetals and others.
In this work, we use a finite element method [46,47] to better present the physical picture of nonlinear rectification via a self-consistent solution of transport equations related to different material parameters. Different harmonics are coupled since the nonlinear equations are solved. Extensive time-domain simulations are performed to uncover the carrier dynamics in the channel following the hydrodynamic and drift-diffusion models. Typically, the non-resonant THz detection of the device with different structural parameters and different doping concentrations has been discussed. Our results not only provide guidance for further device optimization but can also be extended to other novel electronic systems for studying THz nonlinear rectification, which will be very useful for future research into low-energy photon harvesting techniques.

Concept of Nonlinear Rectification in a Quasi-One-Dimensional Analysis Model
Harmonic generation is a general characteristic of nonlinear systems; the process of nonlinear rectification in MOSFETs for THz radiation also exhibits this feature [4][5][6][7][8]40]. Following Dyakonov and Shur's work, the nonlinear properties of plasma waves due to the oscillation of carrier density and velocity play a central role in THz radiation rectification [4][5][6][7][8]. The 2D electron gas will exhibit interesting hydrodynamic behavior [3,4], and the steady current flow in a FET channel can lead to the growth of plasma waves [4]. The AC voltage generated by the incident THz light will act on the channel of the FET, and the nonlinear current-voltage dependence will generate different orders of signal due to the harmonic generation [42][43][44]. Thus, the physical process of THz nonlinear rectification is described as follows.
The surface concentration n in the MOSFET channel is related to the local gate-tochannel voltage swing above the threshold and can be simply determined by the gradual channel approximation [4][5][6][7][8]: where C is the gate-to-channel capacitance per unit area; q is the electron charge and U = U gc (x) − U T ; U gc is the gate-to-channel voltage; and U T is the threshold voltage. The hydrodynamic Euler equation and current continuity equation can be used to understand the dynamic processes of carriers according to the theory of Dyakonov and Shur [4][5][6][7][8], which can be written as follows: where v is the local electron velocity, ∂U/∂x is the longitudinal electric field in the channel, m is the effective electron mass, τ is the momentum relaxation time, n is the electron concentration, and j n is the electron current density (for the n-type channel). Considering Equation (1), the usual continuity equation should be used to calculate Equation (2), which is as follows: These equations coincide with hydrodynamic equations for shallow water. The growth of plasma waves caused by current flow in the FET channel is similar to the shallow water wave [4][5][6][7][8].
Besides, the boundary conditions are as follows: where V gs is the dc gate-to-source voltage swing, U a is the amplitude representing the intensity of the THz radiation, and ω is the angular frequency of the incident THz radiation. U ac = V gs + U a sin(ωt) is the external AC voltage induced between the gate and source by the incoming THz wave, and j = CUν is the electron current per unit width. Then, the solution of Equations (2) and (4) is described in the following form: where − ν, − U are the time-averaged electron velocity and channel potential, respectively, and ν n , U n vary with time with the frequency nω, where ω is the frequency of the input signal. Since Equations (1) and (2) are nonlinear equations, different harmonics are coupled. It is worth noting that the nonlinear properties of plasma waves can provide the possibility of optical detection [4][5][6][7][8]. However, if the input signal, U a , is relatively small, ν 1 , U 1 , are proportional to U a , while − U − U 0 , − ν, and ν 2 , U 2 , are proportional to U a 2 , etc. Because the incident THz electromagnetic waves are small in energy, U a is a small quantity, and U a 2 is an even smaller magnitude. Finally, the embodied higher-order harmonics are dominated by the second-order harmonics.
To the first order in U a , the following equations are obtained: where u 1 = eU 1 /m and s = (eU 0 /m) 1/2 is the plasma wave velocity. Retaining the secondorder time-independent terms, we find the following: Here, − u = e − U/m and the angular brackets denote the time averaging over the period 2π/ω. The boundary conditions for (9)-(12) follow from (5) and (6): For the boundary conditions given by (14) and (15), the integration of (11) and (12) with respect to x yields the following: Furthermore, the evaluation of ∆u from solving the plasma wave dispersion equation, yields the detector response ∆U = m∆u/e, which is the constant source-to-drain voltage induced by the incoming ac signal [4][5][6][7][8].
Depending on the collision-induced damping of the plasma wave oscillations in the channel, the detection mode is composed of resonant and non-resonant modes, which are determined by the excitation frequency ω and the momentum relaxation time τ.
When ωτ 1, the detector works as a resonant detector. The damping of the plasma oscillations is small under the illumination of incident light. The plasma waves reflect between the two boundaries of the channel, causing continuous plasma wave oscillation. The resonance occurs at a specific frequency ω N = ω 0 (2N + 1) leading to a sharp resonance, ω 0 = πs/2L g , s is the plasma wave velocity and L g is the channel length [4][5][6][7][8]24,25].
On the other hand, the detector works in the nonresonant mode due to the overdamping of the plasma wave when ωτ 1. In this case, the DC voltage between the source and drain is a constant. Compared with the case in the non-resonant mode, the generated photoresponse signal could be improved several times in the resonant mode [4][5][6][7][8]. However, several factors, such as channel length and phonons or impurity scattering processes, could limit the signal, which is usually prominent at room temperature. Therefore, we only focus on the THz detector operating in the non-resonant mode at room temperature [4][5][6][7][8].
Theoretically, Equations (1) and (2) indicate the physics of 2D electrons and plasma waves in a simple form. Nonlinear properties of the 2D electron gas lead not only to the rectification of the incoming electromagnetic radiation but also to the emergence of the signal with the second-order and higher harmonics of the incoming radiation [4][5][6][7][8]. It can be seen from Equation (19) that the photoresponse changes only by a factor of three, even when the parameter ωτ increases from zero to a very high value. Even though the underlying physical process may be different, the voltage-tunability of FETs endows them with unique advantages to reach higher performance.

Simulation Model and Equations
Based upon the simplest description of nonlinear rectification in one-dimensional terms, we would like to introduce more intricate situations of nonequilibrium dynamics when considering realistic devices. In this work, Sentaurus Technology Computer-Aided Design (TCAD L-2016.03-SP2) simulation is explored via the hydrodynamic/drift-diffusion transport model with extensive incorporation of different physical quantities, open circuit boundary conditions, etc., which describes the plasma wave in the MOSFET channel in a self-consistent manner [4,[48][49][50]. The drift-diffusion transport model is used when in thermodynamic equilibrium [50]. When further considering the lattice and electron temperatures, we will use the hydrodynamic model included in the software for the simulation [51][52][53]. During the simulation process, the drift-diffusion model, the Mobility model, the Hydrodynamic model and the RSH (Shockley-Read-Hall) model contained in the simulator are fully used. Here, the basic equations used in the Sentaurus Device are described as follows.
Firstly, the software solves the Poisson equation, which is described as follows: where, φ is the electrostatic potential, q is the elementary electronic charge, ε is the permittivity, n is the electron density, and p is the hole density. Secondly, the software solves the continuity equation. In addition, in this device, the charge carriers are electrons, which is described as follows: ∇· → J n = qR net,n + q ∂n ∂t (21) where, R net,n is the electron net recombination rate, → J n is electron current density. Thirdly, the software solves the drift-diffusion transport equation, which is described as follows: where, ϕ n is the Quasi-fermi potential and µ n is the electron mobility. During the simulation, the drift-diffusion transport equations can, eventually, be extended to the hydrodynamic equations. Fourthly, the thermodynamic transport equation differs from drift-diffusion when the lattice temperature equation is solved. However, it is possible to solve the lattice temperature equation even when using the drift-diffusion model. → J n = −nqµ n (∇ϕ n + p n ∇T) (23) where, p n is the absolute thermoelectric power and T is the lattice temperature. Fifthly, when the electron temperature is also considered in the hydrodynamic model in TCAD, the solution equation is the following: → J n = µ n n∇E c + kT n ∇n − nkT n ∇ln γ n + λ n f td n k n ∇T n − 1.5nkT n ∇ln m n (24) in Equation (24), the first term takes into account the contribution due to the spatial variations of electrostatic potential, electron affinity, and the band gap. The remaining terms in Equation (24) take into account the contribution due to the gradient of concentration, the carrier temperature, eventually T e in the channel, gradients, and the spatial variation of the effective masse, m n . for Boltzmann statistics, γ n = λ n = 1. The thermal diffusion constant f td n defaults to zero, which corresponds to the Stratton model in TCAD [51,52]. Theoretically, the carriers absorb photon energy via an interband or intraband transition under radiation, which leads to the average temperature T e rising above the lattice temperature T. This transport model of Equation (24) takes into account the non-equilibrium state between hot carriers and lattice temperature for different channel sizes, which can be used to study the microphysical behavior of hot electrons [51][52][53].
By self-consistently solving the above physical equations, it allows for the description of the microphysical processes at the interior of the device and ensures the accuracy of the simulation results. It will also provide a convenient way to study the nonlinear rectification process in two-dimensional materials, which opens the door for modification of physical and material parameters [29][30][31][32][33][34][35][36][37][38][39].

Simulation Procedure for Finite Element Analysis
Finite Element Analysis (FEA), by using the numerical calculation method, is an efficient tool to solve all kinds of scientific problems [46,47]. Many problems can be solved by using the finite element method, including electromagnetics, heat transfer, fluid mechanics, stress analysis and other problems [54][55][56]. Meanwhile, it is widely used in scientific fields such as aerospace, automotive, electronics, nuclear science and other fields [50,51]. To date, more and more complex problems in the field of science need to be extracted from their physical models and even consider the interaction of various physical fields. At this moment, FEA can be used to combine a variety of complex physical models and solve complex scientific problems [57][58][59].
The finite element method is used by partitioning a continuous entity into multiple subdomains [54]. In this way, the entity is partitioned into a discrete body with no intersection between each subdomain within the body. The discretized subdomains after partitioning are used as differential elements representing the entity model [60,61]. Separate analysis of each differential element is achieved by solving differential equations created for each element and analyzing the solution results of the differential equations [62][63][64].
The physical field equations that satisfy the behavior of carriers in semiconductors [47] can be solved in TCAD, such as voltage and current magnitude, electron-hole movement and recombination [65,66]. The device is divided into numerous subdomains, and the finite element analysis is performed in each subdomain. Each subdomain of the device is interlinked, such as the contact area between the adjacent subdomains, the type of material added, and the doping concentration of the added material [67][68][69]. The specific simulation steps are described as follows.
The schematic of the silicon-based FETs is shown in Figure 1a. The length of drain L d and source L s is 50 nm, the channel length L g is 100-300 nm, and the SiO 2 insulating layer is 10 nm. Figure 1b shows the main parameters used in the simulation. The dc gate voltage V gs and AC voltage U a sin(ωt) are provided between the gate and source (the source terminal is grounded), and V ds is measured with an open circuit boundary condition. Figure 1c shows the finite element division of the mesh structure. Within a specified period of time, the electron concentration of the channel at x = 0 as the gate voltage linearly increases from 0 to 1 V is shown in Figure 1d. During the gate voltage scanning, electron concentration increased by 16 orders of magnitude. Figure 1e,f show the surface electrostatic potential over the range of gate voltage from 0 to 1 V. These figures show that the channel potential continues to rise, and eventually an antipattern layer is formed. Figure 1e shows the cross-section of the channel at y = 3 nm. Figure 1f shows that the surface electrostatic potential rises, extending along the x-axis plane in the channel, as the gate voltage scans. Theoretically, it indicates that the channel is conducting, which can be further used for the internal physical process of the detector. Next, Table 1 lists the important parameters for silicon materials [54].   Figure 2a illustrates the process of finite element simulation and summarizes the main flow of the device simulation. The device is divided into finite grids, and the equations are solved in each grid considering the boundary conditions. The solution results are output after convergence judgment. Figure 2b shows the output current-voltage characteristics of the MOSFET simulated at the gate voltage Vg = 3 V and different gate lengths.   Figure 2a illustrates the process of finite element simulation and summarizes the main flow of the device simulation. The device is divided into finite grids, and the equations are solved in each grid considering the boundary conditions. The solution results are output after convergence judgment. Figure 2b shows the output current-voltage characteristics of the MOSFET simulated at the gate voltage V g = 3 V and different gate lengths. Meanwhile, the simulation results have the same trend as the previously calculated theoretical results, which further validates that our finite element simulation is reliable [70]. Figure 2c shows the transfer current-voltage characteristics of the transistor calculated at the drain voltage V d = 2 V. One can find that the device threshold voltage (V th ) increases as the L g increases. Specifically, when the L g is 300 nm, the V th is 1.23 V, which is significantly larger than that of 100 nm as attributable to a larger parasitic capacitance. Meanwhile, the simulation results have the same trend as the previously calculated theoretical results, which further validates that our finite element simulation is reliable [70]. Figure 2c shows the transfer current-voltage characteristics of the transistor calculated at the drain voltage Vd = 2 V. One can find that the device threshold voltage (Vth) increases as the Lg increases. Specifically, when the Lg is 300 nm, the Vth is 1.23 V, which is significantly larger than that of 100 nm as attributable to a larger parasitic capacitance. Next, detailed quasi-continuous time-domain simulations are performed to analyze THz mixing in the MOSFET channel. Following Dyakonov and Shur's work, the nonlinear properties of the 2D electron gas led to THz rectification and the appearance of the signal at the second and higher harmonics under THz radiation [4][5][6][7][8]. Figure 3a shows the time evolution of Vds at the gate bias of 0.7 V with a channel length of 100 nm. It is worth noting that the drain output is comprised of harmonics and self-mixing DC components. To extract the dc response voltage, we performed a Fourier transform to obtain the frequency domain characteristics of the output signal amplitude (Figure 3b). It can be observed that a fundamental tone of ω= 0.2 THz (the incident frequency) is accompanied by a distinguishable DC component (Figure 3b, in the inset). Moreover, as shown in Figure 3c, the second harmonic (0.4 THz) can be observed.  Next, detailed quasi-continuous time-domain simulations are performed to analyze THz mixing in the MOSFET channel. Following Dyakonov and Shur's work, the nonlinear properties of the 2D electron gas led to THz rectification and the appearance of the signal at the second and higher harmonics under THz radiation [4][5][6][7][8]. Figure 3a shows the time evolution of V ds at the gate bias of 0.7 V with a channel length of 100 nm. It is worth noting that the drain output is comprised of harmonics and self-mixing DC components. To extract the dc response voltage, we performed a Fourier transform to obtain the frequency domain characteristics of the output signal amplitude (Figure 3b). It can be observed that a fundamental tone of ω = 0.2 THz (the incident frequency) is accompanied by a distinguishable DC component (Figure 3b, in the inset). Moreover, as shown in Figure 3c, the second harmonic (0.4 THz) can be observed. In general, the fundamental harmonics (0.2 THz) and the second harmonics (0.4 TH are in along with the DC component (f~0). As a detector, the DC component can be e tracted as the photovoltaic response and is generally used as a signal for THz detection In agreement with Equations (2), (4) and (18), which describe the plasma wave osc lation with an open circuit boundary condition, we can see a significant asymmetric shi in the electron concentration within one oscillation cycle. The asymmetric shift is cause by the interplay between nonlinear electron oscillation and asymmetric boundary cond tions along the channel, leading to a change in potential distribution after time accumul tion. Besides, it can be found that the asymmetric offset becomes more pronounced as th Lg increases, which is probably caused by the loading effect reported before [71][72][73][74][75].  Figure 4d shows the drain output direct current (DC) voltage signal ΔU with th change of gate-to-source static voltage for different channel lengths Lg. With the increas in Lg, ΔU initially increases, and then it can be seen that ΔU reaches its peak at Lg = 20 nm. At channels longer than 200 nm, the signal starts to decrease as a result of the loadin effect from the access region in the MOSFET device, and the efficiency is reduced at th expense of sub-threshold swing.
Besides, in Figure 4e, we plot the optical response of devices with different Lg from 0.1 to 0.9 THz. It can be observed that the response reaches its lowest value at 0.5 THz. A Lg = 100 nm, the response of the device is smaller than in the other two cases. Howeve such a process is inefficient at higher frequencies due to the adverse effects of parasit capacitance and resistance [70]. By taking into account the level of response and reliabilit it can be concluded that the optimal length is at Lg = 200 nm.
As a detector, the linear dynamic response is another important parameter that d scribes the power-dependent properties of targeting weak objects. A good linear relation ship between DC output and normalized input radiation power ( ) obtained by chang ing the power by more than an order of magnitude, is shown in Figure 4f [40]. Indeed, can be inferred from Figure 2 to 4 that the simulated results display very good output an transfer characteristics, and plasma-mixing with high efficiency is given rise at the drain end under THz irradiation and asymmetric open-circuit boundary conditions. These r sults validate the reliability and feasibility of our approach for designing and optimizin such types of devices, according to the literature [53].
Following the results of Figures 2c and 4d, an interesting phenomenon can be iden tified: the Vth is minimal when the Lg is 100 nm, and the saturation curve bends signif cantly upward as Vd increases. In Figure 4d, the response signal does not grow monoto nously after increasing both Vg and Lg and reaches its peak near the threshold gate voltag at Lg = 200 nm. Such a trend cannot be described by simply following Equation (18), from In agreement with Equations (2), (4) and (18), which describe the plasma wave oscillation with an open circuit boundary condition, we can see a significant asymmetric shift in the electron concentration within one oscillation cycle. The asymmetric shift is caused by the interplay between nonlinear electron oscillation and asymmetric boundary conditions along the channel, leading to a change in potential distribution after time accumulation. Besides, it can be found that the asymmetric offset becomes more pronounced as the L g increases, which is probably caused by the loading effect reported before [71][72][73][74][75]. Figure 4d shows the drain output direct current (DC) voltage signal ∆U with the change of gate-to-source static voltage for different channel lengths L g . With the increase in L g , ∆U initially increases, and then it can be seen that ∆U reaches its peak at L g = 200 nm. At channels longer than 200 nm, the signal starts to decrease as a result of the loading effect from the access region in the MOSFET device, and the efficiency is reduced at the expense of sub-threshold swing.
Besides, in Figure 4e, we plot the optical response of devices with different L g from 0.1 to 0.9 THz. It can be observed that the response reaches its lowest value at 0.5 THz. At L g = 100 nm, the response of the device is smaller than in the other two cases. However, such a process is inefficient at higher frequencies due to the adverse effects of parasitic capacitance and resistance [70]. By taking into account the level of response and reliability, it can be concluded that the optimal length is at L g = 200 nm.
As a detector, the linear dynamic response is another important parameter that describes the power-dependent properties of targeting weak objects. A good linear relationship between DC output and normalized input radiation power (U 2 a ) obtained by changing the power by more than an order of magnitude, is shown in Figure 4f [40]. Indeed, it can be inferred from Figures 2-4 that the simulated results display very good output and transfer characteristics, and plasma-mixing with high efficiency is given rise at the drain-end under THz irradiation and asymmetric open-circuit boundary conditions. These results validate the reliability and feasibility of our approach for designing and optimizing such types of devices, according to the literature [53].
Following the results of Figures 2c and 4d, an interesting phenomenon can be identified: the V th is minimal when the L g is 100 nm, and the saturation curve bends significantly upward as V d increases. In Figure 4d, the response signal does not grow monotonously after increasing both V g and L g and reaches its peak near the threshold gate voltage at L g = 200 nm. Such a trend cannot be described by simply following Equation (18), from which only a constant voltage drop can be derived, manifesting the discrepancy of gradual channel approximation from the more realistic case that our numerical modeling captures, which means that the physical processes are more complicated when the device is operated below the threshold voltage V th . Furthermore, our results are also in good agreement with other experimental results for gate voltage-dependent properties, which all show the maximum response near V th [76][77][78].

Impact of Hot Electrons on Optical Detection
In order to understand the role of non-equilibrium hot electronic processes, we have investigated the physical behavior of carriers in MOSFET short-channel structures. When the channel length of the device is less than the thermal relaxation length of the material, the hot electron effect cannot be neglected in the case of detectors [79].
In 2014, Wang et al. first revealed the physical mechanism of photocurrent generation inside the detector, confirming that THz irradiation can generate a current from hot electrons [80,81]. Under the low energy of terahertz irradiation, it forms hot carriers in the channel with an average temperature T e above the lattice temperature T [82,83]. Theoretically, the formation of hot carriers is much faster than the rate of thermal equilibrium with the lattice, and the hot carrier temperature can remain above the lattice temperature for a considerable period of time [79,82,83]. As a result, it creates a temperature gradient ∆T in the channel. During this period, some of the hot carriers in the channel spontaneously diffused from the high-temperature region to the low-temperature region under the drive of the temperature gradient, which generates current in the channel [79,82,83]. For our simulated process, the hot carriers are hot electrons, whose free diffusion leads to hot electron effects that result in a reduced optical response for THz nonlinear rectification. Theoretically, the ability of hot carrier diffusion is expressed by the Seebeck coefficient, and the potential built into the channel under THz radiation is that ∆U = −S∆T [79]. The direction of the resulting current depends on the variation of the material Seebeck system, which is non-monotonic under THz irradiation [80][81][82][83].
The following figure plots the distribution of different electron temperatures during one THz radiation oscillation cycle. In the hydrodynamic model, we have considered the impact of temperature induced by hot electrons on the current. Moreover, other variables are kept constant unless specifically stated. Figure 5a shows that the optical response is somewhat reduced when considering the physical processes in which hot electrons are present. Therefore, we further analyze the internal variation of the detector during a 0.2 THz oscillation cycle. The electron temperature distribution at different times within a THz oscillation cycle (from 0 to 5 ps) is shown in Figure 5c-g. One can find that the highest temperature point of the electron gradually moves from the source to the drain and finally returns to its original state. We further extract the one-dimensional electron temperature profiles along the horizontal direction of the device (Figure 5b), which intuitively illustrates the movement of the electron temperature peak. In this case, a temperature gradient in the channel is formed due to the non-uniform temperature distribution, which results in the directional diffusion of hot electrons from the drain to the source and a reduction of the non-linear terahertz response at the drain output. Following the previous device simulations, our finite-element model can be furth verified by varying the source/drain doping concentrations and examining the hot ele tron effect on the output and transfer characteristic curves of the device [84]. Figure 6a-d display the transfer and output characteristic curves by changing th doping at the source and drain, respectively. At first, we calculate the Vth of the transf curve at different doping concentrations, and Vth is 1.05 V at drain doping of 10 17 cm −3 , V is 1.05 V at drain doping of 10 18 cm −3 , Vth is 1.07 V at drain doping of 10 19 cm −3 , Vth is 1.11 at drain doping of 10 20 cm −3 , and Vth is 1.12 V at drain doping of 10 21 cm −3 . The dopin concentrations in the range of 10 18 to 10 20 cm −3 show an obvious effect on the transferrin curve, while little effect is observed with doping below or beyond this regime. Following the previous device simulations, our finite-element model can be further verified by varying the source/drain doping concentrations and examining the hot electron effect on the output and transfer characteristic curves of the device [84]. Figure 6a-d display the transfer and output characteristic curves by changing the doping at the source and drain, respectively. At first, we calculate the V th of the transfer curve at different doping concentrations, and V th is 1.05 V at drain doping of 10 17 cm −3 , V th is 1.05 V at drain doping of 10 18 cm −3 , V th is 1.07 V at drain doping of 10 19 cm −3 , V th is 1.11 V at drain doping of 10 20 cm −3 , and V th is 1.12 V at drain doping of 10 21 cm −3 . The doping concentrations in the range of 10 18 to 10 20 cm −3 show an obvious effect on the transferring curve, while little effect is observed with doping below or beyond this regime.
Next, we would like to reveal the effect of doping at the source on the device's performance. The V th is growing up by increasing the concentration, e.g., V th is 0.87 V at source doping 10 17 cm −3 , V th is 0.97 V at source doping 10 18 cm −3 , V th is 1.07 V at source doping 10 19 cm −3 , V th is 1.11 V at source doping 10 20 cm −3 , and V th is 1.12 V at source doping 10 21 cm −3 . The output characteristic curves exhibit a saturated trend at intermediate doping concentrations between 10 17 and 10 20 cm −3 . This can reduce the effect of hot electrons on the output characteristics of the device.
Following the above results, the response curves of the device under THz irradiation by changing the doping at the source side are simulated and shown in Figure 7. Under different gate voltage modulations, it can be clearly seen that the device hosts a good response when the doping is between 10 19 and 10 20 cm −3 in Figure 7a. Furthermore, the highest value of the device response shifts toward a larger V g as the doping concentration increases. When the doping concentration is below or above the range from 10 18 to 10 20 cm −3 . Based on this, the optimum source doping concentration can be selected. Next, we would like to reveal the effect of doping at the source on the device's performance. The Vth is growing up by increasing the concentration, e.g., Vth is 0.87 V at source doping 10 17 cm −3 , Vth is 0.97 V at source doping 10 18 cm −3 , Vth is 1.07 V at source doping 10 19 cm −3 , Vth is 1.11 V at source doping 10 20 cm −3 , and Vth is 1.12 V at source doping 10 21 cm −3 . The output characteristic curves exhibit a saturated trend at intermediate doping concentrations between 10 17 and 10 20 cm −3 . This can reduce the effect of hot electrons on the output characteristics of the device.
Following the above results, the response curves of the device under THz irradiation by changing the doping at the source side are simulated and shown in Figure 7. Under different gate voltage modulations, it can be clearly seen that the device hosts a good response when the doping is between 10 19 and 10 20 cm −3 in Figure 7a. Furthermore, the highest value of the device response shifts toward a larger Vg as the doping concentration increases. When the doping concentration is below or above the range from 10 18 to 10 20 cm −3 . Based on this, the optimum source doping concentration can be selected. Furthermore, the optical responses related to different source doping concentrations and different frequencies from 0.1 to 1.0 THz are presented in Figure 7b. It can be found that the optimal doping concentration is between 10 19 and 10 21 cm −3 , and the response remains high with bandwidth exceeding 0.65 THz, which may help in other novel materials Furthermore, the optical responses related to different source doping concentrations and different frequencies from 0.1 to 1.0 THz are presented in Figure 7b. It can be found that the optimal doping concentration is between 10 19 and 10 21 cm −3 , and the response remains high with bandwidth exceeding 0.65 THz, which may help in other novel materials for studying THz nonlinear rectification.

Conclusions
In this work, we have demonstrated the physical mechanism of nonlinear THz rectification by a finite element method by self-consistently solving the transport equation dynamically. Extensive time-domain simulations are performed to show the carrier dynamics in the channel, considering the hydrodynamic and drift-diffusion models. The competitive behavior between nonlinear rectification and hot electrons induced by the photothermoelectric effect is clearly presented in nanoscale detectors. It is found that the hot-electron effect can be reduced by optimizing the source doping concentrations. Our results provide opportunities for device optimization and understanding the nonlinear THz rectification in other electronic systems, enabling low-energy photon harvesting techniques.